Diffusion of muonium and hydrogen in diamond 
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Jump rates of muonium and hydrogen in diamond are calculated by quantum transition-state 
theory, based on the path-integral centroid formalism. This technique allows us to study the influence 
of vibrational mode quantization on the effective free-energy barriers AF for impurity diffusion, 
which are renormalized respect to the zero-temperature classical calculation. For the transition 
from a tetrahedral (T) site to a bond-center (BC) position, AF is larger for hydrogen than for 
muonium, and the opposite happens for the transition BC — > T. The calculated effective barriers 
decrease for rising temperature, except for the muonium transition from T to BC sites. Calculated 
jump rates are in good agreement to available muon spin rotation data. 
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Hydrogen has been studied for many years as an im- 
purity in solids. Due to its low mass, this impurity has 
posed some challenging problems to both experimental- 
ists and theorists [l[ . An important property of hydrogen 
in insulating and semiconducting materials is its ability 
to form complexes and passivate defects [3] , as was found 
in the case of diamond A large amount of informa- 
tion has been obtained by studying muonium (formed by 
a muon fi + and an electron), which behaves as a light 
isotope of hydrogen, with « m p /9. In particular, 
muon spin rotation (/iSR) experiments have provided us 
with invaluable information on the behavior of muonium 
(Mu) in semiconductors 

An interesting topic is the diffusion of Mu and H in 
crystalline materials. This problem turns out to be diffi- 
cult due to the combination of quantum effects with lat- 
tice relaxation around the impurity (polaron effect). In 
fact, zero-point motion influences the excitation energy 
from the ground state to the top of the diffusion bar- 
rier. Also, tunneling of the impurities can be enhanced 
by phonons, and lattice distorsions may depend on the 
isotopic mass of the impurity. 

Two types of isolated Mu centers have been observed 
in elemental semiconductors, and are characterized by 
their isotropic or anisotropic hypcrfine interaction. The 
former is accepted to consist of Mu at an interstitial tetra- 
hedral (T) site, while the latter corresponds to Mu at a 
bond-center (BC) site, midway between two nearest host 
atoms. Muon implantation experiments in diamond, as 
well as various theoretical approaches, have shown that 
this impurity is metastable at the T site, and has its 
lowest energy at or around the BC site, as a result of a 
large lattice relaxation [l], [H, 0] . An important fraction 
of the implanted muons form Mux, and the transition to 
Mubc has been observed by /LtSR Q. This transition T 
— ► BC, and the opposite BC — * T are also expected to be 
important for the diffusion of hydrogen in diamond, as 
derived from the energy barriers calculated with several 
theoretical methods [7|. 



In diamond there is the additional problem that H and 
Mu may exist in different charge states, and moreover a 
change of state may occur in combination with hopping. 
Density functional (DF) theory calculations predict H + 
to be stable off- axis in a buckled bond-centered configu- 
ration, in contrat to the BC site for neutral H 0,0] . Also, 
calculated diffusion barriers were found to change appre- 
ciably with the impurity charge state @, 0] • Here we will 
concentrate on high-resistivity diamond, where the non- 
paramagnetic fraction does not play an important role 
and thus we will study neutral H and Mu. 

In earlier works, hydrogen diffusion in diamond 
has been studied theoretically in the classical (high- 
temperature) limit @, 01- However, quantum effects are 
important for hydrogen-related defects in this material at 
relatively low temperatures Thus, we study in this 
Letter the diffusion of H and Mu in diamond by explic- 
itly considering the impurities as quantum particles. The 
main question to be answered is the dependence of hop- 
ping rates on both temperature and impurity mass. In 
this respect, there is a vast literature about theoretical 
models for quantum diffusion of light particles in solids, 
and metals in particular (loj . Due to the complexity 
of this problem, such computations have been typically 
based on model potentials for the impurity-lattice inter- 
actions. 

Here we calculate the jump rate of hydrogen and 
muonium by quantum transition-state theory Il| , using 
path-integral molecular dynamics (PIMD) simulations. 
We employ a realistic interatomic potential, derived from 
DF theory calculations. This method allows us to cal- 
culate jump rates for this nonlinear many-body prob- 
lem, including lattice relaxation, zero-point motion, and 
phonon-assisted incoherent tunneling. 

In the path-integral formalism of statistical mechan- 
ics, a quantum particle can be represented as a cyclic 
chain of L beads coupled by harmonic springs (L, Trot- 
ter number). This formalism has been employed ear- 
lier to study equilibrium properties of H and Mu in sili- 
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con |l2| and diamond [13j, by using Monte Carlo sim- 
ulations. In this context, there exists a quantum ex- 
tension of classical transition-state theory for calculat- 
ing rate constants of infrequent events [llj. It relates 
the jump rate k to the probability density of the center- 
of-gravity (centroid) of the quantum paths of the jump- 
ing impurity, and particularly to the ratio P c between 
the equilibrium probability of finding the centroid at a 
saddle-point (say r*) and at a stable site (say r ) 
Namely: k = vP c /2l, where I is the distance between 
Tq and r* and v is a weakly temperature-dependent fac- 
tor: v — <f> Ax/((3h). Here, Ax is the width of the 
probability distribution for the jumping impurity with its 
centroid x fixed at r*, f3 = 1/ksT, and 4> is a number of 
order one at low T. P c can be written as exp(— /3AF), 
AF being an effective free-energy barrier, given by the 
reversible work done on the system when the impurity 
centroid x moves along a path from rg to r* : 



AF 



J Tn 



f(x)dx , 



(1) 



where f(x) is the mean force acting on the impurity 
with its centroid fixed on x at temperature T: f(x) = 
— (V x V^(R))x. Here V^(R) is the potential energy, R be- 
ing in our case a 3(N + l)-dimensional vector (N host 
atoms plus one impurity) . The reliability of this method 
to calculate free-energy barriers and jump rates was dis- 



cussed in 1 1 , [1 



We use the Born-Oppenheimer approximation to de- 
fine a potential energy surface V(R) for the nuclear co- 
ordinates. Since true ab-initio potentials require com- 
puter resources that would enormously restrict the size 
of the simulation cell, we have found a compromise by 
using an efficient tight-binding (TB) Hamiltonian, based 
on DF calculations 15j. With this interaction potential 
we found the BC site to be the absolute energy mini- 
mum for H [l3| , and the energy surface is similar to that 
derived from earlier DF and TB calculations @ 

To sample the configuration space we employed the 
PIMD method in the NVT ensemble [H, [ig|. Simula- 
tions were carried out on a 2 x 2 x 2 supercell of the 
diamond face-centered cubic cell with periodic boundary 
conditions, including 64 C atoms and one impurity. To 
assure the right convergence of the path integrals, we 
took a Trotter number L oc 1/T, with L = 20 for H and 
60 for Mu at 300 K. For given impurity mass and tem- 
perature, the mean force f (x) was calculated at 14 points 
along the integration line between ro and r* in Eq. (JlJ. 
For a centroid position x, the nearest and next-nearest 
neighbors of the impurity (up to a total of 13 host atoms) 
are treated quantum-mechanically to obtain f(x), as in 
171 ]. For each point in the integration paths, we gener- 
ated 5000 configurations for system equilibration, and 3 
x 10 4 configurations to calculate ensemble avera ge prop - 
erties. More technical details are given elsewhere 13lll6|. 
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FIG. 1: Effective free-energy barrier for impurity jumps from 
a T site to a neighboring BC site as a function of tempera- 
ture. Open triangles, muonium; filled circles, hydrogen; open 
squares, classical limit. Dashed lines are guides to the eye. 



In Fig. [T] we present the free-energy barrier AF for 
impurity diffusion from a T to a BC site. Data derived 
from line integration of the mean force are shown as a 
function of temperature for hydrogen (solid circles) and 
muonium (triangles). In this plot, one notices first that 
AF is higher for H than for Mu. At low temperature, 
the dependence of AF upon impurity mass is related to 
the change in internal energy F of the defect complex 
along the diffusion path. To compare Mu and H, we note 
that Fmu is always larger than Fh, but the difference 
Fmu — Fh changes from 0.96 eV at a site T to 0.83 eV 
at the transition state, thus giving (AF)h — (AF)m u = 
0.13 ± 0.01 eV. Second, in the case of H, one observes a 
slight decrease in AF for increasing T, contrary to the 
rise in effective free-energy barrier for muonium. At high 
T, one converges to the other, as expected for the classical 
limit. These migration barriers are on the order of that 
obtained earlier from DF-TB calculations in Ref. [S| (0.4± 
0.1 eV). 

For the opposite jump (impurity from BC to T), the ef- 
fective barrier is higher for Mu than for H. This is shown 
in Fig. [2l where symbols correspond to AF derived from 
our PIMD simulations at several temperatures. (Note 
the different vertical scales in Figs. [1] and O) In this 
case, Fmu — Fh changes at low temperature from 0.64 
eV at a BC site to 0.83 eV at the transition state, giv- 
ing: (AF) H - (AF)m u = -0.19 ± 0.01 eV. The low- 
temperature energy barriers obtained here for H and Mu 
are comparable to those found from DF theory (1.6 eV) 
and earlier DF-TB calculations (2.0 ±0.1 eV) [8]. Our 
simulations yield, however, an important decrease in AF 
as temperature is raised. 
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FIG. 2: Free-energy barrier for impurity jumps from a BC 
site to a neighboring tetrahedral T site as a function of tem- 
perature. Triangles, muonium; filled circles, hydrogen; open 
circles, classical limit. Dashed lines are guides to the eye. 



This impurity transition between BC and T sites is 
asymmetric in a double sense. First, it happens between 
a local minimum of the energy surface and the abso- 
lute energy minimum. Second, the lattice relaxations 
involved in both impurity positions are very different. 
At BC, the nearest host atoms relax strongly (~ 0.4 A), 
whereas at a T site the relaxation of the C atoms is much 
weaker (0.08 A). This second asymmetry is relevant for 
the temperature dependence of the free-energy barriers 
shown in Figs. [T] and (2[ For diffusion across static bar- 
riers, it is known that AF increases with temperature 
[l8j |. In our case, the barrier AF depends on lattice 
ralaxation and vibrational modes, and the larger the re- 
laxation, the more inportant is the change in AF with 
temperature. For the transition BC — > T, AF is con- 
trolled by the energy surface around BC, in which the 
lattice relaxation changes appreciably from the BC site 
to the transition point, and thus AF decreases fast for 
rising T. On the contrary, AF for the jump T — > BC 
is controlled by the path between T and the transition 
point, which does not involve important host-atom re- 
laxations, and AF changes slowly with T. In particular 
for increasing T, AF decreases for H but rises for Mu, 
as a consequence of the smaller lattice relaxation for Mu 
(closer to a static barrier). 

In Fig. Owe show the muonium jump rate from a T site 
to a neighboring BC site, as derived from the probability 
P c given by the calculated barriers AF (open squares). 
One observes a certain departure from linearity in the 
Arrhcnius plot, due basically to the change in effective 
barrier as a function of temperature (see Fig. [T]). The 
solid line in Fig. [3] displays the transition rate derived 



FIG. 3: Rate for impurity jumps from a T site to a BC site. 
Open squares represent results derived from PIMD calcula- 
tions. The solid line corresponds to data derived from ^tSR 
measurements |19|. and an open triangle indicates the rate 
measured in [20]. For comparison, we also present PIMD re- 
sults for the jump rate from BC to T (circles). Error bars 
of the simulation data are on the order of the symbol size. 
Dashed lines are guides to the eye. 



by Odermatt et al. [191 ] from /zSR experiments. A tri- 
angle shows the spin relaxation rate of Mut in insulat- 
ing diamond containing vacancies, and measured at room 
temperature [20j. In that work, it was found a rather 
constant relaxation rate at T < 100 K, which could in- 
dicate the appearance of coherent tunneling between T 
sites. For comparison, we also give in Fig. [3] the Mu 
jump rate from BC to T sites, as derived from the free- 
energy barriers shown in Fig. [2] At room temperature, 
and even at T ~ 1000 K, this jump rate is several orders 
of magnitude smaller than that for the transition T — > 
BC. From the jump rate derived from our simulations at 
300 K (fcMu = 3.4 x 10 6 s _1 ), one expects that a muon 
implanted in diamond at a site T will likely diffuse to a 
BC site before decaying (with mean lifetime r M = 2.2 /is). 

We now turn to the diffusion of H in the diamond bulk. 
There are a number of possible migration trajectories, as 
suggested by earlier theoretical works Q. For neutral 
hydrogen, according to the energy surface derived in our 
calculations, one can think of two main diffusion paths 
for hydrogen. The first one (denoted as path I) consists 
of a jump from BC to a neighboring T site, followed by 
a transition to another BC site. This is the process envi- 
sioned in Ref. d from (classical) locally-activated Monte 
Carlo simulations. For both atomic jumps, we obtain the 
effective free-energy barriers shown in Figs. [Hand [2 An 
alternative path (called II) moves H from a BC site to 
a nearest BC site, having a transition state at a point 



4 



with C2v symmetry close to the so-called C site • This 
path is similar to that employed earlier in a calculation 
of the jump rate of H in silicon 11711 . but in diamond the 
energy barrier is much higher [2 11 ] - In fact, for classical 
point atoms at T = we find a barrier of 2.1 ± 0.1 eV, 
close to 1.8 eV obtained from DF calculations Q and 
1.9 eV derived from semi-empirical cluster calculations 
[22j ]. Concerning these energy barriers, it is worth noting 
that seemingly simple atomic jumps can actually involve 
coupled barriers, as clearly indicated in Ref. 23. Thus, 
to obtain the barrier in path II we considered a coupled 
motion of H and the C atom lying between both BC sites. 

Going to the results of our finite-temperature simula- 
tions of H in diamond, we find for path II at 1000 K a 
jump rate from BC to BC sites of 160 s" 1 , much lower 
than that obtained in path I for the transition from BC 
to T (k = 2.5 x 10 5 s _1 ). Taking into account the rate 
for the opposite process T — * BC (k = 4.6 x 10 10 s _1 ), 
and the relative residence time of hydrogen at BC and 
T sites, we obtain at 1000 K a diffusion coefficient along 
path I of D H = 3.1 x 10~ 12 cm 2 s _1 . This indicates 
that H diffusion from BC to BC sites will mainly happen 
via short visits of tetrahedral T sites, or their associated 
attraction basins in configuration space. 

Experimental investigations on hydrogen diffusion in 



diamond have been so far scarce. Popovici et al. [24 1 
studied the diffusion of several species in diamond at 1130 
K, and found that the diffusion coefficients for N, O, and 
H are very similar. They concluded that the diffusion 
process of these impurities is probably affected by crystal 
defects that would trap the diffusing species. Hence the 
value for hydrogen found by these authors [Z?h = (2.4 ± 
0.3) x 10~ 13 cm 2 s _1 ] has to be considered as a lower 
limit for diffusion in a perfect diamond crystal. 

In summary, path-integral molecular dynamics simu- 
lations provide us with a good tool to study quantum 
effects on the jump-rate of muonium and hydrogen in 
diamond. Renormalization of the classical diffusion bar- 
riers due to these effects is appreciable. In particular, we 
have found that the effective free-energy barrier for muo- 
nium can be smaller or larger than that for hydrogen, 
depending on the diffusion process under consideration. 
For muonium, we find jump rates from T to BC sites on 
the order of /iSR data. For hydrogen, the most probable 
diffusion path involves BC and T sites. 

These calculations were performed at the Barcelona 
Supercomputing Center (BSC-CNS) . This work was sup- 
ported by M.E.C. (Spain) through Grant FIS2006-12117- 
C04-03. E.R. Hernandez is thanked for helpful discus- 
sions. 
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